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ABSTRACT 

We present a model of gamma-ray emission from core-collapse supernovae originating from the explosions 
of massive young stars. The fast forward shock of the supernova remnant (SNR) can accelerate particles by 
diffusive shock acceleration (DSA) in a cavern blown by a strong, pre-supernova stellar wind. As a funda- 
mental part of nonlinear DSA, some fraction of the accelerated particles escape the shock and interact with 
a surrounding massive dense shell producing hard photon emission. To calculate this emission, we have de- 
veloped a new Monte Carlo technique for propagating the cosmic rays (CRs) produced by the forward shock 
of the SNR, into the dense, external material. This technique is incorporated in a hydrodynamic model of an 
evolving SNR which includes the nonlinear feedback of CRs on the SNR evolution, the production of escaping 
CRs along with those that remain trapped within the remnant, and the broad-band emission of radiation from 
trapped and escaping CRs. While our combined CR-hydro-escape model is quite general and applies to both 
core collapse and thermonuclear supernovae, the parameters we choose for our discussion here are more typical 
of SNRs from very massive stars whose emission spectra differ somewhat from those produced by lower mass 
progenitors directly interacting with a molecular cloud. 

Subject headings: acceleration of particles, shock waves, ISM: cosmic rays, ISM: supernova remnants, mag- 
netic fields, turbulence 



1. INTRODUCTION 

Many core-collapse supernovae are expected to explode 
within their parent molecular clouds. Because of the influ- 
ence of the surrounding material, the manifestation of the su- 
pernova remnant (SNR) can differ substantially depending on 
the progenitor star type. For a relatively low progenitor mass 
below ~ 12-14 Mq, the stellar wind and photoionizing radia- 
tion are not sufficient to substantially clear out the surround- 
ing cloud and already at a radius of about 6 pc the remnant 
can enter a radiative phase wit h a shock directl y interacting 
with the molecular cloud (e.g.. IChevalierlfl999l) . The radia- 
tive shock with a typical velocity below ~ 150km s _1 can 
acce lerate and compres s CRs and produce non -thermal radia- 
tion dBvkov et alJl2000tlUchivama et al.ll2010l) . Recently, the 
Large Area Telescope on board the Fermi Gamma-Ray Space 
Telescope detected GeV emission from SNRs IC 443, W44, 
and 3C391 kno wn to be directly in t eracting with molecul ar 
clouds (see, e.g jAbdo etaill2010blli[Castro & Slanell2010l) . 

Higher mass young stars with masses above ~ 16 M (of 
BO V type and earlier) are likely to create low-density bub- 
bles and HII regions with radii ~ 10 pc surrounded by a mas- 
sive shell of matter swept up from the molecular cloud by the 
wind and the ionizing radiation of the star over its lifetime. 
In this case, a strong supernova shock propagates for a few 
thousand years in tenuous circumstellar matter with a veloc- 
ity well above 10 3 km s _1 before reaching the dense massive 
shell where it decelerates rapidly. 

Regardless of the SN type, the blast wave of the SNR is 
expected to accelerate ambient material and generate rela- 
tivistic electrons and ions, i.e., cosmic rays (CRs), which 
produce strong non-thermal radiation. A preponderance 
of evidence suggests that the particle acceleration mecha- 
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nism most likel y responsible is diffusive shock accelera- 
tion ( DSA) (e.g..lBlandford & Eichlerf fl987l: iJones & Ellison! 
[l99UlMalkov & Drury||205"li) ~ 

We note that despite the common acceleration mechanism, 
the appearance of the two classes of SNRs we mention can 
differ very substantially. For early progenitor stars, one can 
expect that a sizeable fraction of the 7-ray emission is pro- 
duced by the CR ions that escape the forward shock and in- 
teract with the dense surrounding shell, while for lower mass 
progenitors, the bulk of the non-thermal radiation is likely to 
come from trapped CRs. 

While the CRs produced by the SNR generate non-thermal 
emission across the spectrum from radio to TeV 7-rays, the 7- 
rays are of particular interest because they may be produced 
in proton-proton (or heavier ion) collisions of ultra-relativ- 
istic particles. In fact, there are three populations of shock 
accelerated CRs that are important for producing 7-rays: rel- 
ativistic electrons producing 7-rays through inverse-Compton 
and non-thermal bremsstrahlung; CR ions that remain trapped 
within the forward shock precursor; and CR ions that are ac- 
celerated by the forward shock but escape upstream. These 
three populations are produced simultaneously by DSA but 
they have very different properties and will have very differ- 
ent 7-ray signatures. 3 

As has been known for 
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large fraction of the energy in particles accelerated at strong 
shocks can escape at an upstream boundary. In fact, the 
fraction of all galactic CRs that originate as escaping particles 
is likely to be significant and escaping CRs may even provide 
the bulk of CRs at the "knee" and above. The importance of 

3 A fourth particle population that we don't consider here are secondary 
electro n-positron pairs produced by proton-p roton interactions (see, for ex- 
ample, Gabici, Aharonian & Casanova 2009). These leptons will produce 
inverse-Compton emission and may be important depending on the external 
mass concentration. 
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modeling e scaping CRs was discusse d before the advent of 
DSA (e.g.. lSchwartz & SkilUng|[T978l) and is attracting con- 
sider able attention currently within the DSA paradigm (see, 
e.g.. IPtuskin & Zirakashvilil 120051: ICaprioli. Amato & Blasil 
I2010t lDruryll2010l) . Recently. iReville. Kirk & Duffvl «l2009h 
used a simple iterative scheme to construct stationary numer- 
ical solutions to the coupled kinetic and hydrodynamic DSA 
equations. The stationary solutions with efficient acceleration 
were found when the escape boundary was placed at the 
point where the growth and advection of strongly driven, 
non-resonant waves where in balance. For that particular 
case, they derived the energy dependence of the distribution 
function close to the energy break. As we shall argue below, 
some additional factors, e.g., stochastic Fermi acceleration 
on long-wavelength fluctuations, can affect the spectral shape 
of the escaping particles. 

A number of stationary, nonlinear (NL) models of DSA can 
provide the integrated escaping CR energy flux as a fraction 
of the parameterized overall acceleration efficiency, but no 
model is yet able to determine the spectral shape of escap- 
ing CRs taking into account the self-consistent production of 
magnetic instabilities produced by both the trapped CRs in the 
shock precursor and the escaping CRs. 4 Such a treatment is 
not yet feasible, creating an important problem since the in- 
terpretation of the 7-ray emission from young SNRs depends 
critically on the uncertain spectral shape of both the trapped 
and escaping CRs. Therefore, a suitable parameterization of 
the shape of the escaping CR flux is needed to allow compar- 
isons with 7-ray observations of young SNRs in the hope of 
constraining NL DSA models. 

It is important to note, of course, that SNRs are not sta- 
tionary and the dynamics of an expanding remnant, even in 
the simplest spherical case, adds an additional factor to the is- 
sue of CR escape. As the remnant expands, the precursor re- 
gion beyond the forward shock that is filled with CRs expands 
producing a "dilution" of CR energy density. This effect 
has been studied in detail by Berezhko an d co-worker s (e.g., 
Berezhko. Elshin. & Ksenofontovl Il996alfbh (see also iDruryl 
20101) . In a real shock, the dilution effect is coupled to es- 
cape since the lowering of the CR energy density results in 
less efficient generation of magnetic turbulence and this will 
change the escape of CRs. Both the flow of energy out of 
the shock by escape and the dilution of the CR energy den- 
sity influence the Rankine-Hugoniot conservation relations in 
similar ways. Both act as energy sinks and both result in an in- 
crease in the shock compression and other nonlinear effects. 
In the stationary, plane-shock approximation for DSA used 
here, we ignore dilution and only include the escape of CRs 
at an upstream free escape boundary. It has been shown, how- 
ever, that this plane-shock approximation gives essentially the 
same results for the shock structure as in a spherical, expand- 
ing shock if the specific mode of escape is unimportant (see 

4 In principle, particle-in-cell (PIC) simulations can solve this problem 
exactly. However, it must be noted that PIC simulations cannot yet solve 
the full NL shock problem for SNRs. Most current effo rts with PIC simu- 
lations have been directed tow ard relativistic shocks (e.g.. lSpifkovskvir2008l : 
IRiquelme & Spitk ovskv 2010). Shocks in SNRs are nonrelativistic and non- 
relativistic shocks are harder to simulate than relativistic ones. The NL ac- 
celeration of particles, both electrons and protons, that will produce radi- 
ation spanning radio to 7-rays, requires an extremely large dynamic range 
that no PIC simulation can yet achieve (at least not in three-dimensions), 
and these simulations will not be able to produce results that can be com- 
pared to broadband continuum emission from SNRs for the foreseeable fu- 
ture. Approximate methods, such as we describe here, are needed (see 
Vladimirov, Bvkov & Ellison 2008, for a full discussion). 



iBerezhko & Ellisonll 19991) . The mode of escape becomes im- 
portant, however, if the escaping CR flux is used to calculate 
7-ray emission in material external to the shock. Since we ne- 
glect dilution, our escaping CR fluxes, and the 7-ray emission 
we calculate from them, are over estimated. 

In this paper, we present a new Monte Carlo technique for 
propagating escaping CRs and calculating their 7-ray pro- 
duction via pion-decay, in the circumstellar medium (CSM) 
surrounding the outer SNR blast wave. This treatment of 
escaping CRs is added to our CR - hydro simulation (e.g., 
lEllison. Decourchelle & Ballelll2005t Ellison et alj|2010i and 
references therein) producing a coherent model where a num- 
ber of related elements of the SNR are treated more or less 
self-consistently. The hydrodynamic simulation couples the 
efficient production of CRs to the SNR evolution, including 
the production of escaping CRs as an intrinsic part of the 
DSA process. The escaping CRs are emitted from the for- 
ward shock as the SNR evolves, their propagation is followed 
as they diffuse in the CSM, and the 7-ray emission of these 
CRs is calculated consistently with that from the CRs (protons 
and electrons) that remain trapped within the remnant. While 
a number of important approximations are required, including 
the neglect of precursor CR dilution, this model represents a 
fairly complete and internally self-consistent description of a 
SNR interacting with a non-homogeneous CSM. 

The importance of freshly made CRs interacting with 
their local environment to produce 7-rays has been recog- 
nized for some time and an extensive li terature exists in this 
field. I n a generalization of the model of iGabici & Aharonianl 
(120071) (and previous work, e.g. , lAharonian & Atoyanlll996l) . 
IGabici. Aharonian & Casanova! (120091) calculate the broad- 
band emission, from radio to TeV 7-rays, from CRs pro- 
duced by a SNR interacting with a nearby molecular cloud. 
They emphasize that, depending on the parameters, the 
7-ray emission can exceed other bands by a large fac- 
tor, suggesting that some unidentified TeV sources might 
be associated with clouds illuminat ed by nearby SNRs. 
IGabici. Aharonian & Casanova! (120091) also note the impor- 
tance of the shape of the 7-ray spectrum for identifying GeV- 
TeV sources. 

The model used b y IGabici. Aharonian & Casanova! (120091) 
is based on that of iPtuskin & Zirakashvilil (120051) and in- 
cludes a description of the evolution of the SNR and the 
spectru m of escaping CRs. In IGabici. Aharonian & Casanoval 
d2009l) . the important parameter p max , the maximum cut- 
off momentum for the CR spectrum, is parameterized as 
Pmax(i) oc t~ 5 , where t is the age of the SNR and 8 is taken 
to be oj 2.48 to match the CR data below the knee, as mea- 
sured at Earth. Our parameterization of p max differs consid- 
er ably from this as we discuss below . An important result 
of IPtuskin & Zirakashvilil (120051) (see IBerezhko & Krymskiil 
19881 for an earlier derivation) that is incorporated in the 
Gabici. Aharonian & Casa nova (2009) model and is not mod- 
eled here is that, when integrated over the whole Sedov 
phase, the total CR sp ectrum is a power law of the form 
Fes c oc p~ 4 . Recently. ICasanova et al.l ([20 lOt) have applied 
the IGabici. Aharonian & Casanoval d2009l) model to SNR RX 
J17 13.7-3946 taking into account the details of the ambient 
gas distribution. 

The model presented h ere is similar to that of 
iLee. Kamae & Ellison! (120081) . In both cases, the evolving 
SNR is modeled with a spherically symmetric hydrodynamic 
simulation where the efficient production of CRs via DSA 
is coupled to the remnant dynamics. The main difference 
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is in the treatment of the diffusion of escaping CRs in 
the region beyond the SNR forward shock. The work of 
iLee. Kamae & Ellison! ((2008) uses a "boxel" technique 
whereby, at each time step and each spatial grid in the 3D 
simulation box, particles are exchanged between the adjacent 
boxels according to the particle momentum, location, and 
density gradient. In the model presented here, we use a 
Monte Carlo technique to propagate escaping CRs in the 
region beyond the forward shock. These two methods of 
propagation have distinct advantages and disadvantages, and 
both differ importantly from more analytic models based on 
a direct solution of a diffusion equation. In any case, we 
feel the problem of CRs produced by relatively young SNRs 
interacting with dense, local material is important enough to 
be considered with a variety of complementary techniques. 

Other differences betwe en the boxel model of 
ILee. Kamae & Ellison! d2008[) and our new Monte Carlo 
model are all based on recen t refinements of the CR-hydro 
model (see Ellison e t al.l2010l and references therein) and on 
how we parameterize the escaping CR distribution, described 
below. While these refinements are important for modeling 
specific remnants, they do not substan tially change the results 
given in lLee. Kamae & Ellison! ( 120081) . 

2. MODEL 

The model we present here consists of two main parts. 
The CR-hydro part is used to calculate the evolution 
of a SNR and i s esse n tially the same as that d e scribe d 
in lEllison et"ai~lj2007l) . iPatnaude. Ellison & Slanel (l2009h . 
lEllison et al.1 (1201 Oh and references therein. The evolution 
of the spherically symmetric remnant is coupled to the effi- 
cient production of CRs and the produc tion of thermal and 
non-thermal emission is calculated (see lEllison et al.l 1201 Ol 
for recent work modeling the broad-band emission from SNR 
RX J1713. 7-3946). The diffusive shock acceleration is deter- 
mined in the CR-hydro model u sing the sem i -analytic model 
of Blasi and co-workers (e.g.. | Blasil 120021: lAmato & Blasil 
120051: iBlasi. Gabici"&V annoni 2005|k The injection scheme 
for this model has been discussed in detail in a num ber of pre- 
vious papers (see ICaprioli. Amato & Blasil 120101 for recent 
extensions of the model) but we note that we use a slightly 
different injection method than typically used by Blasi and 
co-workers. Since the diffusion approximation upon which 
the semi-analytic model is based doesn't apply to thermal par- 
ticles, a parameter, r)- m j, must be defined that specifies what 
fraction of thermal particles obtain a superthermal energy and 
are injected into the DSA mechanism. Given this parame- 
ter, the nonlinear DSA mechanism determines the fraction of 
shock ram kinetic energy that goes into superthermal parti- 
cles, i.e., the acceleration efficiency £dsa- The only differ- 
ence in our implementation of this injection model and that 
of Blasi and co-workers is that we specify £ dsa and then set 
77i n j accordingly. Both schemes are approximations since, in 
an evolving SNR, both ifmj and £dsa are likely to be func- 
tions of age. For simplicity, we hold £dsa constant. 

The Blasi et al. model that we use also implicitly assumes 
that the shock is planar and stationary. Apart from the neglect 
of CR dilution, 5 this approximation will be reasonably accu- 
rate as long as the diffusion length of the highest energy CRs 

5 We note that, as for other aspects of DSA, CR dilution will depend impor- 
tantly on the propagation/acceleration model assumed for the highest energy 
CRs. The exact modeling of the highest energy CRs is not yet feasible and we 
parameterize all escape effects with our single parameter / s k defined below. 



10 4 5 rripC, consistent with most models of CR 



is a small fraction of the shock radius. T he sharp X-ray syn- 
chrotron edges ofte n seen in SNRs (e.g.. IWarren eTallEOOl 
lEriksen et al.l 1201 lb implies the presence of amplified mag- 
netic fields which will result in short diffusion lengths. In our 
models here we assume that the diffusion length of protons 
with maximum momentum p max is 1/10 of the shock radius, 
a small enough value to validate the planar approximation yet 
allow p 

production in SNRs 

Accounting for escaping CRs is essential in efficient DSA 
and escaping CRs are implicitly included in Blasi's semi- 
analytic description. However, until now we have not in- 
cluded them in the production of radiation in the remnant 
environment in our CR-hydro model. The neglect of radi- 
ation produced by escaping CRs is justified if the SNR is 
in a uniform CSM with no external density enhancements. 
In this case, the emission from trapped CRs interacting with 
the shocked material is always much greater than that pro- 
duced by escaping CRs in the less dense, unshocked external 
medium (see Model B in Fig. [5j. 

The second and new part of our model is a calculation 
of the escaping CR distribution that emerges from the SNR 
forward shock and the propagation and interaction of these 
escaping CRs in a dense, spherically symmetric shell exter- 
nal to the SNR. Depending on the density of the external 
material, 7-rays produced by the escaping CRs can over- 
whelm those produced by trapped C Rs, as emphasized by 
iGabici. Aharonian & Casanova! (120091) . We note that while 
here we restrict ourselves to spherical symmetry for the ex- 
ternal mass distribution, it is straightforward to generalize the 
Monte Carlo technique to arbitrary mass distributions. 

2.1. Escaping CR Distribution 

As we make clear in describing our parameterized escaping 
CR model, both the fraction of energy in escaping CRs and 
their spectral shape are uncertain. However, while controver- 
sial for some years, the idea that some fraction of the most en- 
ergetic particles in a shock undergoing DSA must escape, re- 
gardless of whether the shock is stationary or not, is now gen- 
eral ly accepted a lthough certain qualifications are still made 
(see lDmrvlBOTol) 



We believe that energetic particle escape is a funda- 
mental and unavoidable part of DSA that must occur in 
all supercritical collisionless shocks regardless of p max 
or time evolution because (1) obser vations and model- 
ing of the Earth bow shock (e.g.. IScholer el at.l 1198 
iMitchell et all 119831: lEllison. Moebius & Paschmannl [199 
support escape, (2) particle escape is an i ntrinsic part of 
many particle-in-cell (PIC) simu lations (e.g.. iGiac alone et"aT] 
U997l:lGiacalone & EHisonll2000l) . and (3) DSA requires self- 
generated turbulence to work over any reasonable dynamic 
range. Since CRs must interact with self-generated turbu- 
lence to be further accelerated, the highest energy CRs far 
upstream in the shock precursor will always lack sufficient 
turbulence to remain nearly isotropic and some fraction will 
escape. These escaping CRs will generate turbulence for the 
next generation of CRs, creating a bootstrap effect. As men- 
tioned above, the dilution of CR energy density that occurs 
in spherical, expanding shocks will be coupled to CR escape 
through the magnetic turbulence generation. 

Given the assumptions and appr oximations of the model, 
the se mi-analytic description of IBlasi. Gabici & Vannonil 
(120051) determines the energy in escaping CRs, Q esc , but does 
not determine the shape of the distribution. While other work 
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does determine the shape (e.g..lviadimirov. Ellison & Bvkov 
2006t IZirakashvili & Ptuskinl 120081: ICapriolT Amato & Blasi 
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2010[) , the shape that results in these models depends impor- 



tantly on arbitrary parameters and the assumptions made for 
the diffusion of the highest energy escaping CRs. 

Since the shapes of the trapped and escaping CR dis- 
tributions, at the highest accelerated energies, are critical 
for modeling both X-ray synchrotron emission and GeV- 
TeV 7-ray emission, we feel it is important to have a flexi- 
ble, i.e., parameterized, model that can be compared to ob- 
servations to provide information on the uncertain plasma 
processes until an adequate theory of self-generated turbu- 
lence in the presence of escapin g particles is developed (see 
iBvkov. Osipov & Ellison! |201 ll for recent work on long- 
wavelength instabilities that may influence the maximum mo- 
mentum CRs can obtain in a given shock). 

As an example of the complexities that may exist, 
the amplified long-wavelength fluctuations discussed in 
IBvkov. Osipov & Ellison! (I201 ll) may result in particle accel- 
eration by the resonant second-order Fermi mechanism. The 
stochastic acceleration rate, t^ 1 , for particles with spatial 
diffusion coefficient n{p) in the shock precursor is t^ 1 oc 
V ph/ K (p)> wnere v ph is the phase velocity. While this rate 
may be below the first-order acceleration rate, it may still be 
high enough to influence the spectra shape at the highest par- 
ticle energies achieved by first-order DSA. The spectral in- 
dex of particles accelerated by the second-order Fermi mech- 
anism depends on the parameter T ar /T^r., w here T esc is the 
escape time (e.g.. iPetrosian & BvkovH2008l) . In the case of 
resonant stochastic particle acceleration by long-wavelength 
fluctuations, T ac /T csc oc M a /[fci r g (p max )], whe re the char- 
acteristic wave number of the C R instability (c.f.. lBellll2004l: 
IBvkov. Osipov & Ellisonll201 lb is 
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Here, j CT is the mean CR current, r g (p max ) is the CR gyrora- 
dius at p max in the magnetic field B, and M a is the forward 
shock Alfvenic Mach number. 

Therefore, for a large enough precursor CR current, j cr , as 
expected for efficient DSA, the parameter T ac /T esc may in- 
fluence the shape of the CR distribution in the spectral break 
region. For instance, if a shock of velocity y s k produces a 
power-law spectrum of accelerated particles up to some max- 
imum momentum and transfers a fraction rj of the shock ram 
pressure to CRs, then r ac /T csc oc f? -1 (c/V^k) with a weak 
dependence on the particle momentum. The smaller T ac /T esc , 
the larger is the second-order Fermi effect and preliminary 
work (A. Bykov, in preparation) suggests that T ac /T csc < 100 
is needed to see a significant modification of the spectral 
shape. While more exact estimates are difficult, we might ex- 
pect rj > 0.5 and V s k > 5000km s _1 to produce a noticeable 
effect. 

When the shock accelerated particles approach p max , they 
begin leaving the upstream region of the shock and the ap- 
proximate power-law distribution of particles that remain in 
the shock turns over in a fashion that will depend on the dif- 
fusion coefficient of the highest energy particles. Whatever 
the plasma processes are for escaping particles, the shape of 
the escaping distribution, F esc (p), is determined by how CRs 
leave the shock and is, therefore, coupled to the trapped dis- 
tribution. Since no current model of self-generated turbulence 
adequately describes the diffusion of escaping CRs, the diffu- 
sion coefficient is generally assumed to be Bohm-like right up 
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FIG. 1. — The curves in the top two panels show the forward shock ra- 
dius, -Rpg, and speed, V[,k> as a function of remnant age. Both ijpg and 
V a lt are determined directly from the hydro code and include the energy loss 
from escaping CRs and all adiabatic effects. The gradual change in slope 
indicates a broad transition between an ejecta-dominated early phase and the 
Sedov phase at later times. The third panel shows the fraction of SN explo- 
sion energy in all CRs along with the fraction going into escaping CRs. In 
the bottom panel, P max from the CR-hydro simulation (soli d curve) is com- 
pared to that used in Gabici. Aharonian & Casanova (2009) (dashed curve). 
While the curves in this figure extend to 10 4 yr, the simulations discussed in 
the remainder of this paper stop at £snr = 1000 yr, well before the SNR 
is f ully in the Sedov phase. We note that the near identity of the CR-hydro 
and Gabici, Aharonian & Casanova (2009) value of p max at 1000 yr is es- 
sentially a coincidence. The parameters used for these results are listed as 
Model B in Table [T] but the quantities displayed here apply to all of our ex- 
amples. 



to p max and independent of position relative to the shock. 

Here, we parameterize the escaping CR phase-space dis- 
tribution, F esc (p), as a modified parabola centered at 
where p max is the maximum momentum CRs would obtain 
if the acceleration cut off sharply when the upstream diffu- 
sion length, ft(p max )/T4k = ^feb* where 14k is the speed 
of the FS and Lfeb is a free escape boundary. In our SNR 
model, the maximum momentum is determined primarily by 
an arbitrary parameter, / s k, which is the fraction of the shock 
radius equal to the diffusion length of protons with momen- 
tum p max , i.e., £ F ebW = fskR s k{t), where i? sk (t) is the 
radius of the FS at time, t. For all of the examples shown 
here, we set / s k = 0.1; a factor small enough to be consistent 
with the planar shock approximation in the Blasi et al. DSA 
calculation. 6 

Our scheme for determining p max gives a very 
different result from the p arame terization used in 
iGabici. Aharonian & Casanova! d2009l) . as indicated in 

6 We note that at early times, setting the acceleration time equal to 
the age of the remnant may give a lower p ma x in which case this 
value is used (see, for e xample, Berezhko, Elshin. & Ksenofontovl 1 19963 : 
Ellison, Decourchelle & Ballet 2005). 
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FIG. 2. — The top two panels show trapped (black curves) and escap- 
ing CR (red curves) proton distributions for simulations where a C ut and N c 
have been varied. The distributions are summed at the end of the simula- 
tion at <snr = 1000 yr and the escaping CR distributions are those leaving 
the FS before any propagation occurs. In the bottom panel we compare our 
para bola fit with ct cu t = 1 and N c = 30 (red dotted curve) to the result 
from Zirakashvili & Ptuskin (2008) (blue dashed curve). Both curves have 
been normalized to the total energy in escaping CRs. 



the bo ttom panel of Fig. [TJ Gab ici. Aharonian & Casanova! 

argue that magnetic field amplification (MFA) may 
contribute to a strong decrease in p max as a function of time 
since it might be expected that MFA is strongest at early 
times, yielding a large magnetic field and a higher p max . 
As the remnant ages, MFA might decrease, producing a 
stronger time dependence than the standard SNR evolution 
would suggest. Magnetic field amplification is not include 
in the examples we show here. We caution, however, that 
nonlinear feedback may reduc e the fu ll effects of MFA (e.g., 
Vladimirov. Bvkov & Ellisonl 120081: ICaprioli et all 120081 
20091) and we feel it is unlikely that a time dependence as 
strong as assumed by iGabici. Aharonian & Casanoval (12009b 
will be obtained. In any case, our main purpose here is to 
introduce a new propagation tool for escaping CRs and not 
be overly concerned with details that are still subject to active 
research. 

An approximate expr ession for the CRs that remain trapped 
within the SNR is (e.g.. lEllison, Decourchelle & Balletll2005l) 



/trap(p) ~ /sA(»exp 



Pn 



(2) 



where /sa ~ (p/Pmax) is the quasi-power law DSA dis- 
tribution obtained by the standard semi-analytic model, a cut 
is an arbitrary parameter that determines the turnover around 
p max , and p max , as mentioned, is determined by the SNR dy- 
namics and / s k- The distribution of escaping CRs, F csc (p), is 
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FIG. 3. — The black and red curves in the top panel show two spherically 
symmetric CSM density profiles. The blue curve in the top panel shows the 
density of the SNR at 4snr = 1000 yr. The two curves in the middle panel 
show the total mass within a particular radius for the CSM profiles at the 
start of the simulation. The dense shell has a mass ofMtot = 10 4 M Q . 
In both cases, CRs escaping from the FS of the SNR propagate in the CSM 
profiles and, at tsNR = 1000 yr, have the color-coordinated density profiles 
shown in the bottom panel. The parameters for the CSM propagation (e.g., 
Eq. [8) are shown in the bottom panel. The density profiles shown in the 
bottom panel are for a particular momentum near the peak of the escaping 
CR distribution. The irregular variations in the escaping CR densities are a 
result of the stochastic nature of the Monte Carlo propagation. 

parameterized by assuming that it is aparabolain log (p 4 F esc ) 
space, that is, 

log [(p'fF^ip)] = 

-a[log(p')-log(l)] 2 +6, (3) 
where p' = p/p ma , x . Initially, we determine b such that 

./trap (Pmax) — -^csc(Pmax) ; (4) 

which yields b = log^" 1 ) = -0.434. 

The width of the parabola, a, is matched to /trap(p) as fol- 
lows. We determine the momentum p c > p max where the 
trapped CR distribution drops by some factor, l/N c , below 
its value at p max , i.e., 

/trap (Pc) 



l/N c 



(5) 



/trap(p 11 ia.> 

Specifying N c uniquely determines p c . We then obtain a by 
setting: 

^*esc(Pc) 



Fesc (pmax) 



that is, 



a(a cat ) = -log[( P ' c )*N c ]/[(logp' c )Y 



(6) 



(7) 



where p' c = p c /pmax and we have written a(a cu t) to em- 
phasize that the width of the escaping distribution depends on 
the cutoff parameter in the trapped CR distribution. The final 
normalization for F csc is set by the total energy in the escap- 
ing distribution, Q C so which is an output of the semi-analytic 
DSA model. 
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In the top two panels of Fig. [2] we show examples where 
a cut is varied between 0.5 and 2 with N c = 3 and 10. 
All other parameters of the CR-hydro model are the same 
for these examples. In all panels, the black curves are the 
CRs that remain trapped in the SNR, /(p), and these distri- 
butions, unlike the escaping CRs, have undergone adiabatic 
losses during the £snr = 1000 yr age of the remnant. 7 The 
parameters, a cut and N c allow a fairly wide range of shapes 
in the critical region around p max , although it is important 
to note that, in our model, for all reasonable values of a cut 
and N c , the escaping CR distribution is expected to be nar- 
row co mpared to the trapped CRs. This di ffers from the 
work of lGabici. Aharonian & Casanova! (120091). as mentioned 
above, and o flOhira. Murase & Yamazakil(l2010l) who assume 
a power-law form for the escaping CR distribution. 

In the bottom panel of Fig. [2] we compare our parameteri- 
zation (red dotted curve) using a r „ t = 1 and N r = 30 to the 
form presented in Zirakashvili & Ptuskinl (120081) (blue dashed 
curve). Other than a sli ght offset of the peak, this cho ice of 
a cu t and N c matches the lZirakashvili & Ptuskinl (|2008f) result 
quite well. We could have obtained an equally good match 
with a different combination of a cut and N c . The quality of 
this match with a cut = 1 leads us to fix N c = 30 and leave 
a cu t as a single free parameter for the coupled shapes of the 
cutoff in the trapped CRs and the escaping distribution. 

2.2. Monte Carlo Model of Cosmic Ray Propagation 

Given the form for the escaping distribution, we propagate 
the escaping CRs using a Monte Carlo technique. 8 As the CR- 
hydro simulation evolves, F osc (p) is calculated for spherical 
shells at time-steps, At, as the forward shock overtakes fresh 
circumstellar material. As the outer-most shell is formed, es- 
caping CRs leave the shell and diffuse into the CSM with a 
momentum and density dependent mean free path given by 

AcSM = A C SM,o('' 9 /' , g,o) Qrg (^CSM/"o)^ /3 " • (8) 

Here, r g = pc/{eB) is the gyroradius, ncsM is the CSM pro- 
ton number density, and a Ig and f3 n are parameters. For scal- 
ing, we use n = lcm" 3 , r g , = 10GeV/(ei?cSM,oX an d 
-Scsm.o = 3/xG. The normalization of the CSM diffusion 
coefficient, Z?csm,o = Acsm,o c/3, can be estimated from 
CR propagation studies (see, for exam ple , iPtuskin et al.l200(H 
iGabici. Aharonian & Casanova! 120091) . For example, with 
-Dcsm.o = 10 27 cm 2 s" 1 , ncsM = 0.01cm" 3 , a rg = 0.5, 
and P n = 1, Acsm ~ 1 pc at 1 GeV, consistent with the fits of 
IPtuskin et all d2006l) . 

We note that the CSM diffusion resulting from Eq. (|8) is 
very different from the diffusion we assume to occur within 
the SNR. For the acceleration process at the FS, we assume 
Bohm diffusion with A ~ r g . Once the trapped CRs have 
been accelerated in the outer shell, these CRs are assumed 
to remain in the shell as it convects and evolves within the 
remnant. In all cases, the CSM scattering is much weaker 
than within the SNR and the escaping CRs quickly fill the 
CSM out to the end of the simulation box. 

7 We note the distinction that the trapped distributions, f(p), in these plots 
determine /sa exactly from the CR-hydro model and the semi-analytic DSA 
calculation, as opposed to the approximate expression for /trap used in Equa- 
tion to fit the modified parabola. 

8 Many of the elements of our Monte Carlo propagation model are sim- 
ilar to that used to m odel nonlinear DSA and are descri bed in detail in 
Jones & Ellison 11991) and Ellison, Barine & Jones 11996) and references 
therein. 



The process continues until £snr is reached during which 
time some number of CR filled shells have been formed 
within the SNR. The escaping CRs fill the CSM region with 
a distribution that depends on Eq. dS) and the properties as- 
signed to the CSM. 

2.3. Circumstellar Medium Properties 

We model the spherically symmetric CSM with a dense 
shell sitting on a low-density, uniform background of density 
n un i. The shell has a maximum density n s h. e U an d an inner 
radius i? s hoii which, for the examples in this paper, is greater 
than the outer radius of the SNR at £snr, that is, the blast 
wave of the SNR has not yet reached the dense shell at the 
end of the simulation. An additional parameter is the total 
mass in the shell, M s h e ii. 

The red curves in the top two panels of Fig. [3] show 
the density and mass distribution for a CSM with n un j = 
0.1cm" 3 , n s heii = 100cm~ 3 , i? she ii ~ 10 pc, and A/ shcU = 
10 4 M Q . Note that the dense shell smoothly rises from n un j = 
0.1 cm" 3 and the rise is centered on i? s hcii- The black curves 
show the CSM with no shell. The total extent of the simula- 
tion box for these examples is ~ 20 pc. Also shown in the top 
panel (blue curve) is the density profile of the SNR at the end 
of the simulation, i.e., at £snr = 1000 yr. The FS, contact 
discontinuity, and RS can be easily discerned from the figure. 

The addition of the escaping CR distribution requires ad- 
ditional parameters and Table Q] gives the parameters for the 
CSM diffusion and the parameters for the external medium. 
As mentioned above, we restrict ourselves to a spherically 
symmetric CSM in this first presentation of the Monte Carlo 
propagation model. 

3. RESULTS 

We use the following environmental parameters for all of 
our examples: (i) the SN explosion energy, E'sn = 10 51 erg; 
(ii) the ejecta mass, M e j = 1.4 M ; (iii) the distance to the 
SNR, c?snr = Ikpc, and; (iv) the ambient magnetic field 
throughout the CSM, B C sm = 3 pG. 

For the diffusive shock acceleration of trapped and escap- 
ing CRs, we fix the following: (i) the fraction of FS radius 
used to determine p max , / B k = 0.1; (ii) the magnetic field am- 
plification factor, £? amp = 1, i.e., no MFA is used; (iii) the 
matching factor defined in equation (0, iV c = 30, and; (iv) 
the DSA efficiency, £dsa = 50%. 

The parameters for DSA and the CSM propagation that are 
varied for our examples are given in Table Q] Again we note 
that we are not attempting a detailed fit to any particular rem- 
nant and that our model is not restricted to the particular val- 
ues for parameters we use here. Any of the environmental or 
DSA parameters can be modified to match a specific object. 

In Fig. [3] we show results for Models A (with a dense ex- 
ternal shell) and B (no external shell), as li sted in Table Q] 
The top two panels were discussed in Section |2~31 In the bot- 
tom panel of Fig. [3] we show the escaping CR densities at 
igNR. = 1000 yr. The escaping densities shown are for a sin- 
gle momentum near the peak of F csc and the parameters as- 
sumed for the CSM propagation are noted in the figure. The 
escaping CRs are emitted from the SNR as it evolves so the 
escaping CRs that were produced earliest have been diffus- 
ing for approximately tsNR = 1000 yr and many have left the 
simulation box. 

For the case where the CSM is uniform (black curves), the 
escaping CRs diffuse outward and uniformly fill the region 
beyond the SNR forward shock with a density that deceases 
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FIG. 4. — The top panel shows the total CR spectra within the forward 
shock at tgNR = 1000 yr (black solid curve: protons; black dotted curve: 
electrons), along with the escaping CRs. The red curves in both panels show 
the escaping CR distribution at the FS, while the blue solid and dashed curves 
show the escaping CRs after diffusing in the CSM profiles shown in Fig. [J] 
The solid blue curves (Model B) are for the constant CSM and the dashed 
blue curves (Model A) are for the dense shell. The two models in this plot 
have a rg = /3„ = 0.5 and Dqsm,o = 1 X 10 27 cm 2 /s. 




Log 10 MeV 

FIG. 5. — Photon spectra at Earth for the examples shown in Fig. [4] As 
expected, the case with a dense external shell (Model A) shines much more 
brightly in 7-rays than the case with a low density, uniform external CSM 
(Model B). The individual components for synchrotron, inverse-Compton, 
bremsstrahlung, and the pion-decay emission from CRs that remain trapped 
in the SNR are indicated. 

uniformly with radius as expected. With the dense shell, the 
escaping CR density drops rapidly as the CRs enter the shell. 
With this M s heii = 10 4 M Q , the shell is about 2 pc thick. Cos- 
mic rays that propagate beyond i? max ~ 20 pc are removed 
from the simulation. 
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FIG. 6. — Escaping CR diffusion in a dense external shell as a function of 
the normalization of Dqsm- In all cases, ct rg = f3 n = 0.5. 

In the top panel of Fig. [4] we show the CR distributions for 
both the CRs that remain trapped in the SNR (black solid 
and dotted curves) and escaping CRs. 9 In all cases, the in- 
tegrated distributions are determined at isNR = 1000 yr. The 
red curve in the top panel is the summed escaping distribu- 
tion as the CRs leave the FS, i.e., before they propagate into 
the external CSM. The solid and dashed blue curves are the 
escaping CRs, at ^snr = 1000 yr, after propagation and we 
remind the reader that our escaping CR fluxes are over esti- 
mates since we don't model dilution which, in fact, occurs 
simultaneously with escape. The bottom panel shows just the 
escaping CRs with an expanded scale. Note that throughout 
this paper we include only escaping protons and ignore escap- 
ing heavier ions and escaping electrons. Trapped electrons are 
considered for inverse-Compton emission. The distributions 
for the escaping CRs after propagation are lower than the dis- 
tribution as CRs leave the FS for two reasons. The first is that 
some CRs escape from the simulation box at i? m ax- The sec- 
ond is that some escaping CRs diffuse back into the SNR and 
these CRs are ignored and not included in the blue distribu- 
tions in Fig. [4] The CRs that remain trapped in the shock are 
summed from the contact discontinuity to the forward shock. 

In Fig. [5] we show the various photon components for the 
models with ct rg = /?„ = 0.5 and Z?csm,o = 1 x 10 27 cm 2 /s. 
The results for the two models are identical except for the 
pion-decay emission from the escaping CRs. As expected, 
escaping CRs interacting with the dense external shell pro- 
duce substantially more emission than those interacting with 
the uniform CSM. In both cases, however, the emission from 
escaping CRs is much more strongly peaked than the pion- 
decay emission from the trapped CR protons. For the trapped 
CRs, the relative intensity of the pion-decay emission and the 
inverse-Compton emission depends on the various parameters 
chosen, most particularly n un ; and K cp , the electron to proton 
ratio at relativistic energies. The fact that our values, n un ; = 
0.1 cm~ 3 and A' cp = 0.01, result in inverse-Compton domi- 

9 In all of the examples in this paper we only calculate CRs accelerated at 
the forward shock and ignore those accelerated by the reverse shock. 
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nating the GeV-TeV emission is not necessarily an indication 
that we believe this will always be the case. The issue is more 
complicated as indicated in a number of recent papers (see, 
for example. iKatz & Waxmartl2008tlMorlino. Amato & Blasil 
I2008L IZirakashvili & Aharonianl l2010t lEllison et al.l 1201 W . 
Regardless of other parameters, the relative importance of the 
7-ray emission from escaping CRs and trapped CRs depends 
mainly on the external density enhancement. 

In Fig. [6] we show the effect of the normalization of the 
CSM diffusion coefficient as escaping CRs diffuse into the 
dense shell. Three effects are noticeable. The first is that the 
escaping CR density drops more rapidly with stronger scatter- 
ing (i.e., smaller -Dcsm,o) as CRs enter the dense shell. The 
second is that the escaping CR density remains larger in the 
region between the FS and the dense shell when scattering is 
strong even though the flux of escaping CRs that leave the FS 
is the same in all three cases. The third effect is that, beyond 
the dense shell, the escaping CR density falls off faster with 
stronger scattering. The CR density remains large within the 
shell (i.e., at radii < 10 pc) for strong scattering because the 
dense shell acts as a valve that slows the flow of CRs out of the 
system. Beyond the dense shell, weak scattering results in a 
more uniform density distribution than strong scattering since 
CRs rapidly fill the available volume when the scattering is 
weak. 

In Fig. [7] we compare the pion-decay emission for the three 
examples given in Fig. [6] all with the same ambient density 
distribution. The CRs trapped in the SNR are the same for 
these cases so the pion-decay emission from the trapped CRs 
(dashed curve) is the same in the three models. Also identi- 
cal for the three cases is the escaping CR flux as it emerges 
from the FS. The sole difference is the scattering strength, 
Z?csm,o> m the CSM and this produces a fairly strong effect 
on the pion-decay emission from the escaping CRs. While 
the emission from escaping CRs shown in Figs. [5] and [7] is 
summed over the entire region from the outer radius of the 
SNR at <snr to i? max ~ 20 pc, when a dense shell is present, 
most of the emission originates in the shell, as expected. 

In Fig. [8] we compare escaping CR distributions for differ- 
ent power-law dependences of the gyroradius, i.e., a rg = 1/3 
(Model E), a rg = 1/2 (Model F), and a rg = 1 (Model G). 
For variety, the models in Fig. [8] along with those in Fig. [9] 
below, use a different set of CSM parameters than the models 
discussed thus far, as shown in Table [JJ For the three exam- 
ples shown, the CSM parameters are identical and the mean 
free paths differ only in the value of a rg ; the normalization 

of the diffusion coefficient I?csm,o = 10xl0 27 cm 2 s _1 and 
(3 n = 0.5 are the same for the three values of a Ig . As the 
dependence on a rg increases, the high momentum CRs are 
able to stream through the CSM quickly and the number that 
remain within the simulation region at ^snr = 1000 yr drops. 
The strong a Yg dependence also results in a flatter radial den- 
sity distribution, as indicated by the blue curve in the bottom 
panel of Fig. [8] The reason for this is that the momentum near 
the peak in the escaping distribution that is used to calculate 
the density profiles is well above 10 GeV so the examples with 
larger a rg have longer mean free paths. 

In Fig. [9] we show a similar plot where we now keep 
a rg = 0.5 and vary the power-law index for the density de- 
pendence of the diffusion coefficient, /?„. When (3 n = and 
there is no density dependence for the diffusion coefficient, 
the presence of the external dense shell produces no effect 
and the red curve in the bottom panel of Fig. |9]falls off uni- 
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FIG. 7. — Gamma-ray emission for the three cases shown in Fig. [6] Since 
only the interaction of escaping CRs with the external CSM is varied, the 
pion-decay emission from the trapped CRs within the SNR is the same for 
the three cases. Referring to Table [T] the black solid curve is Model C, the 
red solid curve is Model A, and the blue solid curve is Model D. 

formly with radius. For stronger density dependences (green 
and blue curves), the density of escaping CRs drops as they 
enter the dense external shell which has a radius i? s hoii — 7pc 
for these models and those shown in Fig. [8] 

In Figs. [lO] and [TTJ we compare pion-decay emission for 
a cut = 1/2 (Model J) and a cut = 2 (Model K). All other 
parameters for these two models are the same as indicated in 
TableUJ As seen in the top panels of these figures, the CR dis- 
tributions vary considerably for these values of a cu t. The pho- 
ton emission (bottom panels), of the trapped (dashed curves) 
and escaping CRs (solid curves) varies less strongly due to the 
fact that the photon emission is naturally spread out partially 
masking the shape of the underlying proton spectrum. One 
clear feature that remains is the low-energy kinematic cutoff 
at a few hundred MeV. Of course, Figs.QJJandQjIwere calcu- 
lated for a particular set of parameters and the relative impor- 
tance of photon emission from escaping CRs versus trapped 
CRs will depend strongly on these parameters. 

4. DISCUSSION AND CONCLUSIONS 

As part of a comprehensive model of an evolving SNR un- 
dergoing efficient CR production, we have presented a Monte 
Carlo technique that describes the diffusion of CRs that es- 
cape from the forward shock of the remnant and propagate 
into a dense, external shell. While a number of calcula- 
tions of escaping CRs and thei r 7-ray production have been 
performed (see, for exam p ie, iLg e^Kamae & Ellison! F2OO8; 
lOhira. Murase & Yamazakil 120101: iDruryl |201Q[ and refer- 
ences therein), there remain many unresolved issues for this 
important problem. Our Monte Carlo method makes differ- 
ent assumptions than analytic calculations based on solving a 
diffusion equation and in some ways is less restrictive, partic- 
ularly if energy losses are included during propagation. 

The important features of our model include: (i) the en- 
ergy content of the escaping CR distribution is determined 
with the shock accelerated CRs that remain trapped within 
the SNR using a planar, stationary, nonlinear model of effi- 
cient diffusive shock accelerat ion that neglects dilution (i.e., 
iBlasi. G abici & Vannonl 120051) : (ii) the acceleration of CRs 
produces changes in the hydrodynamics that modifies the evo- 
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FIG. 8. — Escaping CR distributions, -F C sc, (top panel), and density profiles 
(bottom panel) for Models E, F, and G, as listed in Table [T] The index a Tg 
is varied as shown and n = 0.5 in all cases. For these examples, and those 
shown in Fig. [5] n un ; = 1cm -3 , n s h c ii = 10cm -3 , i? s holl = 7pc, and 
the densities in the bottom panel are for a particular momentum near the peak 
of the escaping CR distribution. The simulation box extends to 12 pc. 
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FIG. 10. — The top panel shows particle spectra for trapped CRs (dashed 
curve) and escaping CRs (solid curve). The bottom panel shows the cor- 
responding pion-decay emission for these distributions along with the sum 
(dotted curve). 
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as shown and cr rg = 0.5 in all cases. 
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lution of the SNR; (iii) the shape of the trapped CR distri- 
bution at the highest energies, which is uncertain due to a 
lack of a well developed theory of turbulence generation for 
anisotropic particles, is parameterized consistently with the 
shape of the escaping CR distribution; (iv) the broad-band 
continuum photon emission from escaping and trapped CRs 
is determined with a single set of environmental and model 
parameters; and (v) although not emphasized or shown in the 
plots here, the thermal X-ray emission is include d consistently 
with the broad-band continuum emission (e.g., lEllison et alj 
1201 Ol and references therein). 

The examples we show indicate the complexity and im- 
portance of including escaping CRs in a consistent fashion 
with CRs that remain trapped within the SNR. The shape of 
the GeV-TeV emission, particularly the low-energy kinematic 
cutoff, is important as one of the main ways of determin- 
ing whether this emission is pion-decay or inverse-Compton. 
If other features are discernable, they may provide clues to 
the importance of the escaping CRs and external density en- 
hancements. We note that all of the spectra shown here are 
integrated over the region between the contact discontinuity 
and the forward shock and are not line-of-sight projections. 
It should be clear from Fig. [6] that line-of-sight projections 
might show additional strong effects as escaping CRs interact 
with nearby dense material. Line-of-sight projections will be 
included in future work. 

An important parameter that we haven't varied here is the 
efficiency of DSA. In all of our examples we set £dsa = 
50%, i.e., 50% of the forward shock ram kinetic energy flux 
goes into CRs (trapped and escaping) at any instant. In fit- 
ting an actual SNR, £usa is a parameter that may or may not 
be constrained by the observations. A considerable amount 
of work has led to the conclusion that £dsa ~ 50% is a 
likely figure for young SNRs but this efficiency will defi- 
nitely vary between remnants, may vary during the remnant 
lifetime, and may even var y at different locations in a sin- 
gle S NR (see, for example, IVolk. Berezhko & Ksenofontovl 
120031) . We note that Monte Carlo shock simula- 
tions that include M FA and have parameters t ypical 
of young SNRs (e.g., IVladimirov. Ellison & Bykovl 120061: 



IVladimirov. Bvkov & Ellisonll2008l) . show total acceleration 
efficiencies £dsa > 50% with a sizable fraction of total shock 
ram kinetic energy (> 30%) placed in escaping CRs. 

Since we set £dsa = 50% for all of our examples, and the 
other SNR parameters that determine what fraction of explo- 
sion energy ends up in CRs are kept constant, the third panel 
in Fig.Q]gives the results for all of our models. After 1000 yr, 
~ 30% of the supernova explosion energy has gone into all 
CRs with — 10% going into escaping CRs. At 10,000 yr, 
~ 50% has gone into all CRs with ~ 20% going into escaping 
CRs. 

In this initial presentation of our Monte Carlo technique, 
we have exploded the supernova in a uniform CSM with an 
external, spherically symmetric shell of dense material. This 
simple scenario shows how important escaping CRs can be 
for modeling non-thermal emission of young SNRs. It is not 
meant to match any particular object. The Monte Carlo prop- 
agation part of the CR-hydro model can be easily general- 
ized to include asymmetric external mass distributions, such 
as those expected when remnants interact with a dense molec- 
ular cloud (e.g., SNR RX J1713.7-3946). Future work will 
also model 7-rays produced when escaping CRs interact with 
the complex structure of a dense surrounding shell as expected 
from a progenitor stellar wind. 
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TABLE 1 

Parameters for DS A and circumstellar medium diffusion. 



Model 




Ctcut 


-DCSM,0 


ACSM,0 


Org 


0n 




"shell 


A/shcll 


-Rshcll 








[cm 2 s~ ] 


[pc] 






[cm -3 ] 


[cm -3 ] 


[M ] 


[pc] 


A 


lxicr 2 




1 x 10 27 


3.3X1CT 2 


0.5 


0.5 


0.1 


100 


10 4 


10 


B 


lxicr 2 




1 x 10 27 


3.3xl0~ 2 


0.5 


0.5 


0.1 








C 


lxicr 2 




1 x 10 26 


3.3xHT 3 


0.5 


0.5 


0.1 


100 


10 4 


10 


D 


lxlO" 2 




1 x 10 28 


3.3X10" 1 


0.5 


0.5 


0.1 


100 


10 4 


10 


E 


lxlO" 2 




1 x 10 27 


3.3X1CT 2 


1/3 


0.5 




10 


10 3 


7 


F 


lxlO -2 




1 x 10 27 


3.3xl0~ 2 


0.5 


0.5 




10 


10 3 


7 


G 


lxicr 2 




lxlO 27 


3.3xl0- 2 


1 


0.5 




10 


10 3 


7 


H 


lxicr 2 




1 x 10 27 


3.3xl0~ 2 


0.5 







10 


10 3 


7 


I 


lxicr 2 




1 x 10 27 


3.3xKT 2 


0.5 


1 




10 


10 3 


7 


J 


lxicr 4 


0.5 


1 x 10 27 


3.3X1CT 2 


0.5 


0.5 




10 


10 3 


7 


K 


lxicr 4 


2 


lxlO 27 


3.3xl0~ 2 


0.5 


0.5 




10 


10 3 


7 



